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ABSTRACT 


A higher order shear deformation theory is used to 
carry out the stress and stability analysis of laminated 
composite plates. A displacement based finite element 
formulation using a isoparametric, 9-mocl<s,ot. loc^ro- 
-r^g'iQn element is employed. A finite element buckling 
formulation for the higher order theory is developed and 
presented. Separate formulations are done for symmetric 
and unsymmetric laminates to achieve computational effi- 
ciency. A FORTRAN program is developed to efficiently 
implement the formulation. Numerous examples drawn from 
available literature are solved to emphasize the reliabilit 
of the formulation and program. The results are presented 
in the form of tables and graphs. 



CHAPTER 1 


INTRODUCTION AND LITERATURE REVIEW’ 

1»1 Introductory Remarks 

The fibre reinforced composites are being increasingly 
used as load bearing structural members in aeronautical, 
automobile, marine, civil and other fields of engineering. 
This is due to the superior mechanical properties possessed 
by these materials. A composite material offers a signifi- 
cant weight reduction in a structure in view of its high 
strength to weight and high stiffness to weight ratios. 

Also, in fibrous composites the mechanical properties can 
be varied by suitably varying the orientation of the fibres. 
Other advantages include improved fatigue life, wear resi- 
stance, corrosion resistance, temperature dependent behaviou 
etc. 

A fibrous composite material essentially consists of 
high strength fibres embedded in a matrix. The most common 
types of fibres that are in use include glass, carbon, 
graphite, boron, silicon carbide, kevlar, etc. A few of the 
matrix materials being used are epoxy phenolics, polyesters, 
etc. 
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A lamina is a thin sheet with fibres oriented in a 
certain direction forming the main load bearing members. 

The matrix is characterized by a low modulus and high 
elongation and it transfers stress to the fibres through 
bonding at the fibre-matrix interface. It keeps the fibres 
in position, provides flexibility and protects the fibres 
from the environment. A lamina can be basically considered 
to be homogeneous, two-dimensional and orthotropic. 

A laminate consists of a number of laminae bonded 
together. The anisotropic behaviour of a laminate can be 
tailored through the variation of laminae stocking sequence 
and fibre orientation. This feature' gives additional degree 
of flexibility to the designer while simultaneously posing 
unique challenges to the analyst. The extensional modulus 
along the fibres is usually very large relative to the 
extensional moduli! in the lateral directions and the trans- 
verse shear moduli!, due to the differences in elastic pro- 
perties between fibre filaments and matrix materials. This 
is a marked departure from conventional isotropic materials. 
Consequently, the relative importance of physical effects is 
influenced by the directianal nature of properties and theii 
relative magnitude. 

Due to the relatively low transverse shear modulii, 
the shear deformations are much significant in fibre rein- 
forced composites. Further the delamination mode of failure 
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is initiated due to interlaminar stresses* The inplane 
stresses govern the strength requirements* The buckling 
phenomenon becomes important for aircraft panels made of 
composite materials, as the buckling of the skin affects th^ 
aerodynamic characteristics giving rise to inadmissible drac 
at high speeds. Therefore, an accurate prediction of the 
structural response with regards to strength and stability 
becomes imperative, 

1.2 Literature Review 

The classical laminate theory, which is an extension oi 
the classical plate theory to laminated plates, ignores the 
transverse stress components. The first complete laminated 
anisotropic plate theory is attributed to Reissner and 
Stavsky Ll]* The classical laminate theory is adequate onl'^ 
for very thin plates. 

An adequate theory for laminates must account for 
accurate distribution of transverse shear stresses. Many 
theories which account for the transverse shear and normal 
str4sses are available in literature, prominent among them 
being that of Mindlin [2] and Reissner [3], A generalizatic 
of the first order shear deformation theory for homogeneous 
isotropic plates to laminated anisotropic plates is due to 
Yang et al. [4j» Whitney and Pagano [5j extended the YNS 
theory [4] to laminates consisting of arbitrary number of 
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bonded anisotropic layers. In this theory f the normals to 
the midplane before deformation remain straight but not 
necessarily normal to the midplane after deformation. Con- 
sequently, a correction to transverse stiffnesses is require 
The shear correction coefficient depends on constituent ply 
properties, ply lay up, fibre orientation and the particulai 
application. Reddy 16 } used a C penalty plate bending 
finite element for the solution of the governing equations 
of YNS theory. The displacements (u, v, w) and midplane 
rotations were independently interpolated and 

coupled via the shear energy terms which behave as penalty 
terms that restrict the shear rotations to zero at thin 
plate limits. Subsequently, Reddy {,1 \ presented a compari- 
son of the finite element solutions using the penalty method 
and exact closed form solutions, for thick anisotropic 
rectangular plates. The effects of reduced integration, 
mesh size and element type on the accuracy of the penalty 
finite element were investigated. He concluded that the 
9-noded Lagrangian element is not affected by integration 
scheme while the 8-noded element is profoundly affected. 

The consideration of symmetry in the analysis of unsymmetri- 
cally laminated anisotropic plate was studied by Reddy [S] . 
The equilibrium equations for a lamina were used by Engblora 
and Ochoa [9] to arrive at a more accurate prediction of the 
transverse shear stresses. Equations were written in terms 
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of differences in stresses across a lamina* cnalonous to th< 
equilibrium equations. Thus a set of simultaneous equationi 
were obtained considering all the laminae and solved. 

The deflection analysis of hybrid laminated rectangula] 
and skew composite plates was conducted by Iyengar and 
Umaretiya [lO]. They used the Galerkin rechnique to obtain 
numerical results for Kevlar/epoxy and Boron/epoxy laminate5 
It was observed that for a given aspect ratio, an increase 
in skew angle increased the rigidity of the skew plate. The'i 
also concluded that for a specified deflection, the hybrid 
laminates are lighter. Haas and Lee [llj used a nine-noded 
degenerate solid shell element based on modified Hellinger- 
Reissner principle. Functions were assumed independently 
for inplane and transverse shear strains. Stress smoothing 
was adopted by Owen and Li [l2j to overcome the discontinuil 
of transverse shear stresses. 

Further developments in laminate analysis have been 
higher-order theories. ’Higher order’ denotes a variation oi 
displacement that is higher than linear through the laminate 
thickness. A discussion of these theories appears in Ref. 12 
Lo et al. [l2,13j and Reissner [14] presented theories for 
plates based on assumed higher order displacement field. The 
displacement field proposed by Lo et al. incorporate inplane 
modes of deformation in the displacement field of Reissner. 
A higher order theory which not only accounts for transverse 
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shear deformation but also satisfied the zero shear stress 
conditions on the top and bottom faces of the plate was 
suggested by Reddy [15]. Phan and Reddy [16] presented a 
finite element formulation of the higher order theory in 
which inplane displacements are expanded as cubic functions 
of the thickness coordinate while the transverse deflection 
is kept independent of the thickness coordinate. The con- 
tributions and effects of transverse normal stress and strain 
are neglected. Putcha and Reddy [17] developed a mixed shear 
finite element based on the higher order theory. The element 
consists of eleven degrees of freedom per node, i.e., three 
displacements, two rotations and six moment resultants. 

The transverse shear stresses and inplane displacements 
are considered to be continuous between layers in Ren*s 
theory [18], He obtained closed form solution for anti- 
symmetric laminates. Krishnamurthy [l9] assumed a displace- 
ment field which leads to parabolic inplane stresses but 
neglects transverse normal stress. The condition of zero 
transverse shear on top and bottom free surfaces was satis- 
fied. Later [20] he used cubic inplane and parabolic trans- 
verse displacements to obtain closed form solution for a 
symmetric four-layered infinite strip. He observes that 
’ statically equivalent’ estimates of interlaminar stresses 
display the required interface continuity and agree very 
closely with exact solutions. VijJaykumar and Krishnamurthy 
[ 21 ] presented an ' iterative modelling’ in which displacement 


i 
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field is successively generated by utilising the displace- 
ment field of immediately previous lower order theory. This 
implements the higher order theories in a hierarchial fashior 
Librescu and Khdeir [22] employed the state space concept to 
obtain solutions for governing equations of Reddy's type 
higher order theory. The effect of various boundary condi- 
tions vjas investigated. Kant [23j based his higher order 
shear deformation theory on the Taylor's series expansion 
of generalized displacements. He used a new approach for 
satisfying zero transverse shear conditions on top and 
bottom bounding planes. A variety of models based on vary- 
ing assumptions for the displacement fields have been deve- 
loped and tested for accuracy. A displacement based C° 
finite element formulation was presented. 

The finite element formulation of the plate buckling 
problem is available in standard texts, such as that of 
Zienkiewicz [24] and Cook [25] and various papers [26,27, 

28]. The buckling of multilayer plates has been dealt with 
by Lun^gren and Salama [29] using the finite element method. 
Buckling analysis of composite plates has been discussed 
in Iyengar's text [30j. Jones [31] obtained closed form 
solutions for the buckling and vibration of unsymmetric 
cross ply laminated plates. A method for incorporating the 
curvature terms in the nonlinear strains used in developing 
the initial stress matrix is given by Whitney [32]. He 
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concludes that curvature terms have little effect on the 
critical buckling loads for the symmetric laminates that 
were investigated. 

The vibration and stability and symmetrically laminated 
composite plates were studied by Dawe and Craig [33] using 
the finite strip and Rayleigh Ritz method. Owen and Li [34] 
used the heterosis element for the vibration and stability 
analysis of laminated plates. Recently, Gajbir Singh and 
Sadasiva Rao [35] used a 8-noded element with five degrees 
of freedom per node for the stability analysis of thick 
angle ply composite plates. They use all the components of 
the nonlinear strains and the initial stress matrix is obtain 
ed by carrying out the stress analysis. The effects of fibre 
orientation, material properties, layering and boundary con- 
ditions on buckling load are studied. 

Phan and Reddy [36] used the Reddy’s higher order theory 
to obtain closed form solutions for the stability of isotro- 
pic, orthotropic and laminated plates, A mixed element 
based on the same theory is developed by Putcha and Reddy 
[3^]. Doong [38] adopts a higher order theory for the vibra- 
tion and stability of an initially stressed laminated plate. 
He compares the results with Reddy's [36] closed form solu- 
tions. Khdeir [39] exactly solves the higher order governing 
equations developed by Reddy, by making use of the generali- 
zed Levy’s type solution in conjunction with the state space 
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concept. He obtains uniaxial critical buckling loads for 
symmetric cross ply laminates for varying boundary condi- 
tions, Subsequently, [40j he used his higher order theory ' 
analyse the free vibration and buckling of unsymmetric cros: 
ply laminates. Recently, Ren [41j obtained results for 
vibration and buckling analysis of simply supported, unsy- 
mmetric, cross ply and angle ply laminates. His theory 
allows for parabolic distribution of transverse shear stres; 
which are continuous across the interface. 

A review of the different approaches used for modellinc 
multilayered composite plates is recently published by Noor 
and Burton [42j. The discussion focusses on the approaches 
for developing two-dimensional shear deformation theories. 

1,3 Objective and Scope of the Present ISIork 

The literature review reveals the following facts 
prompting and influencing the course of the present work : 

* A very few displacement based finite element adoptions 
of the higher order theory are available in the literature. 
This necessitates the assessment of the reliability and 
accuracy of such a formulation for varied physical condi- 
tions and material parameters. This becomes all the more 
imperative due to the simplicity and versatality of the 
finite element method for use in practical problems. 
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■* A comparative study of the numerical results 

obtained from two or more higher order theories in order to 
establish their relative predictive capabilities has never 
been taken up. 

* The adaptability of a formulation for laminated compo- 
sites to the analysis of sandwich plates has been rarely 
tested. And except for Kant [23], none seems to have 
approached the sandwich plate problem through the higher 
order theory. Hence, this area remains largely unexplored 
and untested. For example, the consideration of shear 
rigidity of faces needs to be thoroughly examined. 

* Formulation of a higher order theory for the buckling 
of laminated composites using the finite element method is 
attempted in only one or two papers [37]. Further, the 
displacement based formulation for the same seems untouched 

In the light of the above observations, an attempt has 
been made in the present work to analyze the bending and 
buckling response of laminated composite plates. The higher 
order theory proposed by Tarun Kant [23j forms the basis for 
the present study. The displacement based finite element 
method using a isoparametric, lagrangian, quadrilateral 
element is employed. A consistent formulation for the buck- 
ling analysis is devised and presented. The results are 
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obtained for a variety of examples chosen from available 
literature which cater for the variation of a number of 
parameters. Transverse shear stresses are evaluated using 
constitutive relations as well as equilibrium equations. 

The application to the stress analysis of sandwich plates is 
extensively investigated. The results are compared with 
those from first order shear deformation theory and classi- 
cal plate theory in order to determine the range of applica- 
bility and necessity of each. The present results are also 
compared with the results from other higher order theories, 
particularly closed form solutions. 

Formulations are presented separately for symmetric ans 
unsymmetric laminates in view of the computational efficienc 
achieved thereby. A program is developed in FORTRAN for the 
computer implementation of the formulations. 


CHAPTER 2 


THEORETICAL FORilULATION 

2.1 Preliminary Remarks 

This chapter presents a detailed formulation of a higher 
order shear deformation theory for the bending and buckling 
analysis of laminated composite plates. The displacement 
field is chosen in a series form. The governing equations 
are developed entirely in terms of the displacement at the 
midplane. Strains and stresses are obtained from displace- 
ments using strain-displacement and constitutive relations. 

The development of the theoretical formulations for symmetric 
and L'Miyrmetric laminates are presented side by side. 

2.2 Assumptions made in the Analysis 

(a) Displacements and strains are small and linear elasticity 
theory applies. 

(b) The transverse normal stress and strain are negligible. 

(c) The points on the normal to the reference surface before 
deformation need not remain on a straight line after deforma- 


tion 
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The assumption (c) is the basis for any higher order 
theory* It implies a nonlinear variation of the inplane 
displacements through the plate thickness* 

2*3 Definition of the Displacement Field 


The displacement field chosen is such as to approximate 
the 3-"dimensional elasticity problem by a 2--dimensional plate 
problem. This involves the expansion of the displacement 
components (u,v,w) at any point in the plate in a Taylor's 
series in terms of the coordinate normal to the reference 
surface, *Z'. The inplane displacements are expanded as 
cubic function of thickness coordinate while the transverse 
deflection is assumed to be constant throughout the thickness 
This is in accordance with the assumptions (b) and (c). The 
displacement field is defined as ; 


u(x,y,z) = u(x,yp) + z(|f)o + |r 

i <1 

v(x,y,z) = v(x,y,o) + ^ ^ 


+ 1_ ,3 , A. 




(^) + 1_ 


z^ (^) 

^ 3^0 

dz 


w(x,y,z) = w(x,y,o) , (2.1) 

where x,y and z represent the coordinate/axes of the plate. 


The expansions of the inplane displacements (u,v) acco- 
mmodate the warping of the transverse cross-section. Certain 
terms in the expansions of the displacements given by 


14 


Eqns* (2el) correspond to the membrane behaviour while 
others represent the flexure behaviour* The two contri- 
butions can be grouped and written as follows : 


u = 

V = 

w = 


membrane part 


2 * 

U + 2 U 
O O 


2 * 

+ z 

o o 


flexure part 

3 ♦ 

+ Z © + 2 

+ z © + z^ ©3 

y y 

Wq (2«2a) 


where, u_ = u(x,y,o) j v. = v(x,y,o) and w. = w(x,y,o) are 

o o o 

the displacements at the midplane, ©^^ = ^ ®y ~ ^3z^o 

are rotations of the normals to the midplane about y and x 

axes respectively. The parameters (u., v_) and (9„, ©„) 

u o X y 

are the higher order terms in the Taylor's series defined 
at midplane. They can be identified as the higher order 
contributions to the inplane displacements and rotations 
respectively. 

For symmetric laminate, the midplane of the plate re- 
presents the neutral surface during deformation. Thus, when 
only transverse loading is considered, the membrane parts 
can be dropped from Eqns. (2.2a) to yield a simplified versio 
of the displacement field for such class of laminates, viz., 

3 * 

u = z © + z © 

X X 

v = z ©y + z^ ©y 

w = Wq (2.2b) 

■to 
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The geometry of a laminate plate with assumed positive 
set of coordinates and physical midplane displacements is 
shown in Figure 2.1. The lamina reference axes are also 
defined. 

2.4 Stress-Strain Relations for an Orthotropic Lamina 

Each lamina of a composite structure is treated as 
homogeneous and orthotropic. A laminate is formed by bonding 
together a number of such laminae. 

Noting the assumption that w is constant across the 
thickness of the lamina and that normal stress "the 

stress-strain relationship is obtained. Thus, the stress- 
strain relationship with reference to the principal material 
axes for an orthotropic lamina can be written as. 



In general, it may be required that each lamina has its 


prin 
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respect to the laminate reference system of axes. This 
necessitates the transformation of the above constitutive 
relations from the lamina axes (1,2,3) to the laminate 
reference axes (x,y,z). This transformation is done using 
the relation. 


= t”^ c [t“^]^ 


(2.4) 


where the transformation matrix T is defined as, 


SC 

0 

0 


c 

-sc 

0 

0 


-2SC 

2SC 
2 2 


0 

0 


0 

0 

0 

c 

s 


0 

0 

0 

-s 

c 


(2.5) 


in which C = cosa and S = sina , 


The triple product defined by Eqn. (2.4), now relates 
the stresses and strains for an orthotropic lamina with 
reference to the laminate reference system of axes, given as. 




^11 ^12 ^13 ^ ^ 



%2 ^23 ^ ° 

,cr 

xy 

r = 

Q 33 0 0 

<7 

yz 


%4 %5 



Sym 





or, in concise form, a = Q€ 


y 

i- 

Y 


'yz 


xz 


( 2 . 6 ) 
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The coefficients of the Q matrix are given 


^11 ^ ■*■ ^22 


Ci2 (S^ + C^) + (Cj^1“C22*"4C33) S2 


(Cii~Ci 2~2C33) S + (C, o-Coo+SC^-,) C 


'12 ""22 ^^^3 3 ' 


Cii + C22 + (2Cj^2 


(Cii~Ci2“2C33) S C + (Ct o-Coo+2C-3Cf) S C 


12 ""22^^""33' 


(Cj^2”^‘^12‘''^22”^^33^ ■*■ ^^33^^"^ 


^44 ^ C55 S 


'55 '^44 


C..) S C 


%5 " *^44 ® ^ ^55 


^ij ^ %’i » i,j = 1 to 5 


(2.7) 


where a stands for the orientation of the lamina principal 


axes with respect to the laminate reference axes, as indica- 


ted in Figure 2.1, 


2. 5 Strain-Displacement Relationships 


The definition of strains as per the theory of elasticity 


is adopted in obtaining the strain-displacement relations, 


ed b 


n to 
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for symmetric laminates 



_ du 

s= 



dv 


2 Ky 

«y 

= ^ 



du 

+ 

dv 

^ xy 

= ■57 

^ = 


dv 


dw 

^yz 

SS 

+ 

3y = 

^xz 

du 


dw 

= 




3 * 

r K 

X 


3 


3 ^ 

z K z K 
xy xy 


2 ^ ^ 
: 0 

y 


2 ^ ^ 

- K 


(2,8a) 


where, 


d© d@ d© d© 

(ir V ic \ sz _-X 4, 

y* ^ 6x * by * 3y dx 


(K , K , K ) = (- 53 - — , ^ 


^ * 

d© d© 

X . Y \ 

37 “ 


f 0 f 0 0 "" 

X* ^y* '^x ^ y 


dw dw ^ ^ 

5 = <®x + ST’ ®y * 3 y^- 3®x > 3ey ) 


(2.8b) 


In concise form, the flexural and transverse shear 
strains for any layer can be separately written using 
matrix notations as, _ ^ 


= z < K. 


w 


o ^ 

zK+z'^K 


(2.9a) 
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= < 

‘Vyz’' 

b 1 

!». 4 

y 

, + z^ 

0 """ 
Y j 

». = 0+z^0* 


X 

J 


(Z5 

X 


0 * 

X 



t2.9b) 


for unsymmetric laminates : 


= 3^ = «x. £„ » + z- K 


3 * 

X 


£ =^ = £ +2K +Z^£ * +Z^K* 

^ ^Yn Y Yn y 


^ dv _ ^ 

y xy = ■5y + 33? - ^xy. 


+ z Kxy + 6^y* + z3 K^y 


^ + ^ = 0 + 21|J +z2 0* 


yz uz Oy 


Y =^ 4 .^ = 0 ^-ziIj 4-z'^0 

' XZ OZ Ox X ^X X 


2 ^ * 


(2.10a) 


Some of the above zompcnsnts are already defined in Eqns. 
(2.8b). The other components are defined now as : 


(e,. ^ ^ 

o » y^ , ^xy^ 


\/0 O 0,0\ 

' ^ 6 X * 3y * dy 3x 


and (ijj^, fy) = (2u^*, 2vJ) 


(2.10b) 


The expressions for strain in any layer given by Eqns. 
(2.10a) can be represented in concise matrix form as, 
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= 25 4- 2i1j + (2.11b) 


2« 6 Laminate Constitutive Relations 

The strain energy of the plate, which is basically 
defined in terms of physical stress and strain components 
can be conveniently written in terms of stress resultants 
and strains. The stress resultants so used are nothing but 
the integral of the physical stress components through the 
thickness of the plate. These stress resultants for the 
differential element of a laminate are expressed in terms 
of midplane plate curvatures and shear rotations, which are 
a function of the assumed, displacement field. 

The strain energy *U’ for the plate with volume *V' and 
surface area *A’ can be written as 
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U = T cr dv (2.12) 

V 

Writing in terms of the individual stress components, 

-f£a + Y o 4-Y cy +Y cr Idv 
2 ^ X X y y ' xy xy 'xz xz 'yz yz-* 

(2.13) 

Substituting the expressions for the strain components 
given by (2.8a) and (2.10a), and conducting explicit inte- 
gration through the laminate thickness, the strain energy is 
obtained in terras of stress resultants, i.e., 

U “ -T ^ Of dA (2.14) 

A 

where e and a are defined as, 
for symmetric laminates : 

i” = (K , K , K , K K K *, 0 , 0 , gJ 0 *) 
fc \ X* ‘ y* xy* X * Y xy * y’ x* y * x. ' 

= (K, K^, 0, 0*) (2.15a) 

5- = (M^, My, M^y, m/. My", M^", Sy, S^, Sy", s/) 

= (M, m", S, s") 


(2,15b) 
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for unsymmetric laminates : 


£ 


* 


(^v t » *= 


5cy 


X V 

O 0 


xYo f Kyf ^ 


K 0 , , t|j , f , 0 *, ^Z^ *) 

(Cq, e/, K, k"", 0, 1{J, 0*) 


(2.16a) 


(Nj,» Ny» ^xy’ * ^y * ^xy' * ^x* ^y’ ^xy* 

* # -3^ ^ 

M ,M ,M ,S,S,T,T,S ,S ) 
x» y» xy* y x y* x* y’ x 

(N, N*, M, M*, S, T, S*) 


(2.16b) 


where, the components of the stress resultant vector a are 
defined by 


''i 

{N,. N„. N„„) =■ Z 


x^ Y 


L =1 
NL h 


(M^, My, M^y) = E 

L=1 

t * NL 


"L+1 

■I' <°x> °y’ 
'’l 


L-H 


X +1 

r 

‘l 


^K’ N^y’')= / (\> °y. 


(M^, M„, U )= Z 


X ^ y ' xy 


NL “L+1 


L=1 


.L _3 


f (‘^x* '^y* °^xy^ ^ 

^L 



23 


tSx- 


NL 

= £ 

L=1 

’’l+i 

i 

o ) 



NL 

= £ 
L=1 

’’l+1 

‘'l 

•3 ,) 
yz'^ 

/ * 


NL 

= £ 

L=1 

'’l+1 

'’l 

<J ) 
yz* 


L 

L 

L 


dz 


zdz 

o 

z dz 


(2.17) 


The stresses are now written in terms of the genera- 
lized strain components given by (2.15a) and (2.16a). 

This is accomplished by substituting strain expressions 
given by (2.9) and (2.11) in the expressions (2.6) for the 
stresses. Thus, the stresses in the Lth layer are given 
by the fol lowing relations : 


for symmetric laminates : 


f ^ 

‘’x 

L 


— — 

L 


r 1 

Kx 


— — - 

L 

Kx 

“y 

r 


Q 





Q 

z 

* V 

S f 

xy 






^y 




* 

^vJ 


fcr 1 
yz 

< 

L 

> 

'^44 *^45 

L 


+ 

^^44 *345 

CM 

N 


\z 


_^45 


L<^xj 


i>45 %5_ 




(2.18) 



24 


for unsymmetric laminates : 



These expressions for stresses, when substituted in 
expressions ( 2 . 17 ) yield the laminate constitute relations, 
which relate the stress resultants to the generalized strains 


as follows : 
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wh6ir©yHj^ — ■j' *” ^ ~ l»2j3y o«»a T 

The relations (2.21a) and (2.21b) can be represented 
together in a concise matrix form as follows : 



(2.22a) 

or, a = D e (2.22b) 

Similarly, for the particular case of the symmetric 
laminates, the relations (2.21) reduce to the form as given 
below, where the coupling terms vanish. 



(2.23a) 
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Sy 


^44^1 *^45^1 

'^44*^3 

^45^3 


f 

>- 



%5^1 

Q 45 H 3 

%5^3 


^x 

is ** 
y 

r ~ 

Sym, 

«44H5 

Q 45 H 5 

4 , 

«'y 

* 

Sx 






fx* 





— 


— 


(2,23b) 


2,7 Interlaminar Stresses Using Equilibrium Equations 

The stress-strain constitutive relations of the foim 
(2,18) and (2.19) are generally used for the evaluation of 
interlaminar shear-stresses present theory 

leads to realistic parabolic transverse shear stress distri- 
bution through the laminate thickness. But, when the con- 
stitutive relations (2.18) and (2.19) are used, they lead 
to the discontinuity of interlaminar stresses at the inter- 
face of two adjacent lamina with non-identical Q matrices. 
This violates the equilibrium conditions. In the present 
section, an alternative approach is adopted for the determina- 
tion of o and d . The equilibrium equations of elasticity 
for each layer are used to derive expressions for the inter- 
laminar shear stresses. The expressions for inplane stresses 
from constitutive relations are utilised. 


We have, for each layer, 




4" 


(2«24a) 




d a 


XV 


da da 


5x oy 


ll ^ lll I tM^T i «i» M y SZ 

Ov oz 


In an alternate form, 


xz 


z=h 


L+1 




Z I dz 

i=l h . 


da ^ 

y 

dy 


yz 


z=h 


L+1 


L 

£ 

i=l 


h.+l 

i 

/ ( 




da. ^ 
Tx 


+ -T^) dz 


(2.24b) 


(2.25a) 


(2.25b) 


The expressions for the inplane stresses given by (2.18) 
and (2.19) are substituted in the above and integration is 
performed through the thickness of the layer. This yields 
the following expressions : 

for symmetric laminates : 


xz 


z=h 


- £ [Qj_i'"(H 


L+1 


i=l 


d^© 


d^e d^© * 

^ + H. -^) 


‘ 2 x 2 

dx 


dx^ 


>.2^ * 
d”^© 


^12 ^^2 33 ^ ^ ^4 ^ 


d^© 6^© d^© * _ 

+ Qi3(2H2 33^ + H2 t 2H4 dx<3y ^4 .2 ^ 


^2^ ^ 
d © 


dx^ 
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d^e. 


dV 


-.9 * 


v2 * 


*^33 ^^2 ^2 ^4 ^4 


(2.26a) 


yz 


L i d^e 

z=h|^+l 1=1 ^ 




x2 * 

d^©. 


+ QiaCHg 


d^e 


d^© 


52 © * 


v2 * 

d^© 


+ QljCHg ^ + 2 H 2 33^4 + H^ + 2H^ 3333^) 


d2© 


•v2 ■*• 

d © 


- «i3<H2 " «4 j;;r-> 


o 

d © 


d^©^ 


v2 ^ 

d'^©. 


^\* 

+ Q33(H2 ^x5y '*' ^2 6xdy ^^4 6x<^y ^4 


^)3 


(2.26b) 


For unsymmetric laminates ; 

L 


xz 


z=h 


L +1 


d^u d2© 

• [Clil(Hi ^ H2 ^ + H3 ^ 


i2 * 

o 


+ H 




d^© 


dx^ 

%2 * 

® ®x ^ 
4 

dx 

a2 * 

0 V, 


a2©* 


+ dxdy '*' ^2 dx^y % dxdy ^4 dx^y ^ 


d^u 


d2© 


a2 * 

o uo 


. o # 
d^9 

X 
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a2 * 
0 V 


d2e * 


^ H, ^ + H2 ^ ^ H3 ^ 


+ Qi.(H 


33-^1 5y2 




i2 * 

d'^v 

+ H-3 , 4. 1^ 

3 dy-^ 


. 2 * 
d^9 . 


4 dy2 


d^u d^© d^9 * 

■" ^ ^ ^2 ^ ■" ^3 ^ ^4 


U V w ^ 

^2 ^i3y + % 


i.2 * 

d V 


.2 * 

d^© 


4 dxoy 


L . d'^u 

.5, Ks (Hi 3757 * «2 3^ ^ H3 

a2e * , 

, T, X \ 


‘4 d^y 


a2 * 

0 u 


(2,27a) 


, a2v„ aV aV* 

+ Q22 (Hi ^ + Hj ^ + H3 + H4 ^^2 


a2 * 

0 n 

o 


>2 * 
d*^© 

X 


'%3 (Hi 5^2 H2 ^^2 ■*■ H3 g^2 + H4 ^^2 

d\ d^e 

5x6y ^2 dxdy ^3 6xSy ^4 'Sx(^y" 




s 2 ^ 

> d'^u 

r^ + H« ■ ■ - A —*— + H, 


^2.^ * 
0 © 


+ Q33(Hj_ + H 2 + H3 


a2 * 

d u_ 


* d2e * 

O *L, O' 

V ”4 dx6y 


d\ d^e ^^®v* 


^4 :k 2 

dx 


(2.27b) 



CHAPTER 3 


FINITE ELEMENT DISCRETIZATION 


3«1 Preliminary Remarks 

The finite element displacement formulation is used 
for obtaining the solutions for the fundamental equations 
derived in the second chapter. The continuum is separated 
by imaginary lines into a number of finite elements. A set 
of functions is chosen to define the state of displacements 
over each element in terms of its nodal displacements. The 
total potential of the system is then written in terms of 
the prescribed displacement field and the principle of 
minimum potential energy is applied to yield the equilibrium 
equations . 

3.2 Element Shape Functions 

The domain can be discretized into ’NEL’ number of 
elements and the total potential (n) of the structure be 
written as a summation of the total potential of individual 
elements (is®), i.e.. 


n = 


NEL 

E 

e=l 


It 


e 


(3.1) 
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0 

In the displacement based formulation, % is written 
as a function of unknown displacement variables defined by 
assumed displacement models. For the displacement models 
given by Eqns. (2.2a) and (2.21s) respectively for unsymme- 
trie and symmetric cases, the vector of unknown displacement 
variables can be written as. 

Symmetric j h = (wQ,0j^,©y,©^*,©y*)^ (3.2) 

Unsymmetric : h = 

(3.3) 

The variation of these over the element is defined in 
terms of normalized coordinates ^ and rj, by isoparametric 
shape (interpolating) functions, such that, 

m 

6 = £ N. 6. (3.4) 

i=l ^ ^ 

where, NN denotes the number of nodes per element, is the 
shape function associated with node i, and is the genera' 
lized displacement vector corresponding to the ith node. 

In the present formulation all the components of the genera- 
lized displacement vector of an element h , are defined by 
identical^ shape functions. The nine-noded lagrangian qua- 
dratic plate element is adopted which accounts for parabolic 
variation of the displacements over the element. The shape 
functions are given by, 
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for corner nodes ; 


N® = •! (1 ^ ^i) (1 + f (3.5a) 

for inidside nodes t 

^ i f ^ I T)Ti^(l+Tir)^)(l- (3.5b) 

for the central node 


N 9 = (1 ~ ^^)(1 ” (3.5c) 

By definition, the shape functions for an isoparametric 
element are used to define the coordinates at any point on 
the reference x~y plane of the plate, viz.. 


NN 

S 

i=l 


%x.. 


m 

z 

i=l 


^i^i 


(3.6) 


3.3 Element Stiffness Matrix 

The stiffness matrix is obtained by expressing the 
internal strain energy for an element in terms of the 
unknown nodal displacements. The strain energy is contri- 
buted by membrane, flexure and shear strains, the membrane 
strains being absent for the particular case of symmetric 
laminates. 
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The membrane strains , bending curvature and shear 
strains are written in a matrix form in terms of the nodal 
displacements using the Equations (2.8b) and (2.10b), as 
given below : 


'®o' 


”k " 


“0 " 


^ ^ 1 

^ m 

K* 

— 

II 



lu 1 
0"" 



(3.7) 


* I 

For symmetric laminates G^, G^ and Ip do not exist. 

All the vectors and matrices L^, and L in relation 

b’ s m 

(3.7) are described below ; 


for symmetric laminates : 

o — (w ,0 ,G ,6 ,0 ) 

^ X* y’ X ’ y ' 

(K% K ) = ,Ky ,K^y ) 


( 3 .8a ) 
(3.8b) 



= (0^, 

0 » 0 

0y") 


0 

d/dx 

0 

0 

0 

0 

0 

6/dy 

0 

0 

0 

d/dy 

d/dx 

0 

0 

0 

0 

0 

d/dx 

0 

0 

0 

0 

0 

d/dy 

0 

0 

0 

d/dy 

d/dx 




(3,8c) 


(3.8d) 
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Ls = 


d/dj 


(3.8e} 


for unsymmetric laminates 


^ X 

*^o *®x ’®y ^ 


(3.9a) 


(®o ’®o ^ ^ ~ ^®x ’®v ’®v ®xv ^ 


"O ^ o 


(3,9b) 


(k^k" ) 


(K^,Ky,K^,K/,Ky'',Kj^y'') 


(3.9c) 


(Qf^,l|j^,0* ) = (0jj,0y,t|J^,llJy»(Z5/»0y'') 


(3.9ci) 


d 6 
5y ^ 


L„ = 0 0 

m 


0 0 


0 0 


"Sy 3x 


(3.9e) 



0 

0 

0 d/dx 

0 

0 


0 

0 


0 

0 

0 

0 0 

d/dy 

0 


0 

0 


0 

0 

0 

0 d/dy d/dx 

0 


0 

0 


0 

0 

0 

0 0 

0 

0 


0 

d/dx 

0 

0 

0 

0 0 

0 

0 


0 

0 

d/dy 

0 

0 

0 0 

0 

0 


0 

d/dy d/dx 

0 

0 

d/ dx 

1 

0 

0 

0 

0 

o“ 


0 

0 

d/dy 

0 

1 

0 

0 

0 

0 
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0 

0 

0 

0 

2 

0 

0 

0 
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0 

0 

0 

0 

2 

0 

0 


0 

0 

0 

0 

0 

0 

0 

3 

0 


0 

0 

0 

0 

0 

0 

0 

0 

3 



(3.9g) 


The generalized displacement vector is defined in 
terms of its nodal values in relation (3.4) using inter- 
polating functions. Using this relation, along with the 
expressions for strains in matrix form given by relation 
(3.7), the generalized strain vector can be expressed 
in terms of nodal displacement values as. 



m m 

2 6. = E Bm. = Bm d 

i=l ^ ^ i=l ^ ^ 
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K 


K 


m 


m 


> =Lb b = L^ S N. bi =1 b. =Bb d 


i=l 


i=l 


0 


, m m 

^ u b = £ N. b. = £ b. = d 

’ S S.,11 .,s. 1 s 

1=1 1=1 1 


0 


(3.10a) 


where, B = N. , B^. = N. , B^^ = 

1 11 


(3.10b) 


and 




B 


^!NJ 


®b = f^b^ ^b2 


B. ] 
“^NN 




Bg j 
NN 


(3.10c) 


and 

= (b, , 


,T _ T 1 T 


1 » ^^2 


f # # o < 




(3.10d) 


The internal strain energy of an element can be 
obtained by multiplying the stress resultants with appropri- 
ate strains and integrating them over the area of the 
element, i.e.. 
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U® = I / [(sj, sf ) -i 


+ (0*^, M dA 

«• 


N 

_N* 

r3 


r + (K^, ) 


T IM 




M' 


(3.11) 


The stress resultants can be substituted by midsurface 
strain vectors using Equations (2.22aJ. Thus, the expression 
(3.11) changes to, 


U 


^ i t T f \ r\ I T ^ ri^ 


, / [(6‘, S* ) D„ 
2 ^ o o m 


<1 




in 




+ (K*^, K^^) 

T 1 

^0 

-i- (k"^, K*”^) Djj ^ 

-K^-l 


i?oJ 


Lk*J 


+ (0"^, 0* ) Dg M dA 

f 


'oil 


(3.12) 


Substituting the strain vectors by the nodal displace- 
ment vector, using Equations (3.10a), 


u® = J / [d^'CE^ D^B„)d + d^CB^ D^Bj,)d + d^CBj B^)d 

A 

-s- d'^(B^ 

or, 

U® = I d^ [k]® d 


(3,14) 



[k]® is identified as the stiffness matrix for an element e. 
For symmetric laminates, the quantities in the above 
expression related to coupling do not exist. 

3.4 Geometric Stiffness Matrix 


For the stability problem, the geometric stiffness 
matrix is derived by considering the work done by inplane 
stresses cf^ and while acting through the second 

order contributions from lateral displacement w to strains 

®y ®xy* 

These strain contributions are : 



- i 
~ 2 



f dw\2 


dw dw 

33 mMmmm mraum) 


(3.16) 


The work done by inplane stresses on an element while 


acting through the above strains is : 
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(3,17) 


After integrating through the thickness, the expression 
for work done can be written in terms of stress resultants 
in a convenient form as, 

dw/dxl 

dw/ayj J 

It is well known that a ’consistent* geometric stiffness 

matrix yields more accurate buckling loads and mode shapes 

30 

which are proven upper bounds to the exact solution. In 
this approach, the displacement functions as used for deri- 
ving, the elastic stiffness matrix are used to build the 
geometric stiffness matrix. 

The derivatives of w can be written in a convenient 
matrix form as. 



dw 



g 


(3,19) 


6 is already defined in (3.8a) and (3.9a) respectively for 
symmetric and unsymmetric laminates. The matrix for the 
two cases is written below : 
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for symmetric laminates : 

d/dx 000 




d/dy 0 


0 


0 


0 


(3.20) 


for unsymmetric laminates 




0 0 d/dx 0 0 0 0 0 o' 


0 0 d/dy 0 0 0 0 0 0 


(3.21) 


The derivatives dw/dx and dw/dy can be expressed dx 
in terras of nodal displacements using relation (3.4) as. 




dw 


dw 

Zy 


L„ 6 

g 


m NN 

= Z N4b4 = Z 64 = d 

i=l 


^ -1=1 ^ ^ i.=l ^ ^ 


( 3 « 22 ) 


w 


here, B_ = [b 


„ B„ .... ] and d is as 

gi 92 


defined in (3,10d). 

Substituting (3.22) in (3.18), we get, 


W- = I I [d" bJ 


^xy| 

L^xy 


Bg d] dA 


(3,23) 


or 


W® = I d"^ [Gj® d 


(3.24) 
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where, 

G® = / IB^ 

A ^ 

G®is known as the geometric stiffness matrix corresponding 
to an initial state of membrane stresses, defined by N , N 

j 

and N , . 
xy 

3.5 ESnement Load Vectoi 


N N 
X xy 


N N 
xy y 


dA 

9 


(3.25) 


The external loads considered for the bending problem 
may be independent or combination of the load cases that 
follow : 

(a) concentrated nodal loads (P ), acting in the direction 
of the coixer/jccn Jinc; nodal degree of freedom. 

(b) uniformly distributed load (P^) over the top or bottom 
faces of the element acting in z-direction. 

The total external work done by these forces on an ele- 
ment can be expressed as : 


Wp = d^ [P^ + / P^ dA] 

Pi 

or 

Wp = d^ P® 

P® is the element load vector and P^ 

c 


(3.26) 


(3.27) 


Cl C 2 


TjT 

'NN 


43 


where P is the vector of concentrated forces/couples at 
node i in the direction of nodal degrees of freedom. 

vector is defined as, 



where, P^^ is given by, 
for symmetric laminates : 


Pq = (Pq, 0, 0, 0, 

'^i 

for unsymmetric laminates : 

P^^ = (0, 0, Pq, 0, 0, 0, 0, 0, 0)"^ (3.28) 

3.6 Derivation of the Equilibrium Equations 

The total potential of an element can be written as : 


= U® - 


(3.29) 


Using Equation (3.1), the total potential of the whole 
structure can be written in general as follows : 


NEL 

ts = 2 

e=l 


U® - W 


(3o30) 
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0 

U is the strain energy of the element given by the 
relation (3®14), W® is the total work done by the external 
forces » 

For the bending problem, the total external work done 
by applied forces is denoted as W® given by the relation 
(3,27). Using the expressions (3.14) and (3,27) for strain 
energy and external work done respectively, the total 
potential for the bending case is, 

^ = 2 t(i d^ [K]® d - d^ P®] (3.31) 

e=l ^ 


minimization with respect to d leads to the following 
equilibrium equation system, 


NEL 

Z 

e=l 


iK]® 


d = 


NEL 

2 

e=l 


(3.32) 


which can be solved for the displacement vector'd'. 

For the buckling problem, the total external work done 
by a set of inplane loads is denoted as W®, given by the 
relation (3,24). Using the expressions (3,14) and (3.24) 
for the strain energy and work done, the total potential for 
the buckling case is, 

NEL 

IS = 2 [(| [k]® d) - (| d"^ [g]® d)J (3.33) 

e=l 
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For convenience, N , N and N may be expressed in 

X y xy 

terms of a common load *N' , such that 


= a N 
Ny = p N 

= Y N (3.34) 


Substituting in (3.25) and assuming the initial state 
of stresses to be uniform throughout the plane of the 
element, the geometric stiffness matrix can be written as 

r« yi 


[GJ® = N S IbI 
A ^ 




Bg] dA 


(3.35) 


= N LG’ ]® 


The applied stresses are assumed to increase in propor- 
tion to a load factor »f’ upto the level when buckling 
33 

occurs. Therefore, at the point of buckling. 


[Q]® = fN [Q’]® 


(3.36) 


writing, = fN, 


(3.37) 


[G]® =N,^[G’3® 


(3.38) 


At the point of instability, the total potential has a 
stationary value. Hence, using (3.38) in (3.33) and minimi- 
sing % with respect to d, we get. 
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[[K]® ~ [G* j®] d = 0 


(3.39) 


which is a typical eiO'.iir.^s.lue problem. On solving, the 
lowest magnitude eigenvalue gives the critical buckling 
load. The eigenvectors represent the buckled shape. 


3.7 Evaluation of the Stiffness Matrices and Load Vectors 

The integrals involved in the expressions for [kJ®, 
[G’J® and P® are evaluated using the Gauss-Quadrature. To 
achieve economy in computation, the element elastic stiff- 
ness matrix and the geometric stiffness matrix are parti- 
tioned in both the row and column directions into blocks. 
Each block corresponds to the degrees of freedom at a node. 
The expression for the stiffness matrix for a typical block 
ij can be written as, 


1 1 


K 


3-J 


/ / D B. iJl d§ dp 

-1 -1 ^ J 


or 


K? . 


m n 


2 2 Ip jJl B". D B. 

f=l g=l ^ y ^ J 


(3.40) 


Similarly, 

1 1 T 

G!® = / / bI 

“^3 -1 -1 Vi 


a Y 

Lt PJ 


Bg^ 1 j 1 df dp 


^ NG NG , T 

or, an =2 2 jjj B ' 

f=l g=l ^ ^ ^i 


a Y 
Y ^ 


B. 


(3.41) 



47 


in which and Ig are the weighting coefficients, ‘NG* is 
the number of numerical quadrature points in each of the 
directions § and t|* 1j| represents the determinant of the 

standard Jacobian matrix that relates the shape function 
derivatives in the natural coordinates ^ and tj with those 
of X and y coordinates. The subscripts i and j refer to 
the nodes that are related to block i j . 

The matrices and Bj are given in general as. 



01 ® 





m . 

3 

Bb i 
°3 

II 

m 


and 


< 


®s 

1 





1 

, i 

J 





j 


(3.42) 


where the membrane terms are non-existent for symmetric 
laminates , 


The computation of K. . and Gl . is economised by 

2.3 

explicit multiplication of the triple products. Also, due 
to symmetry of the stiffness matrix, the block K. . lying on 

Awl 

one side of the main diagonal are only formed. Each block 
in Gl . contains only one one-zero term which only is requirec 
to be formed. 

The evaluation of the load vector P® defined in 
Equation (3.27) is also done using numerical integration 
by expressing as follows : 



using Gauss quadrature. 


p® = p + z S ij| ? (3.43) 

f=l g=l f 9 o 

The exact integration of the integrals in the stiff- 
ness matrices requires 3x3 Gauss point integration rule. 
But, at the thin plate limit, this may lead to overstiff 
results. This is because, the shear energy terms in the 
elastic stiffness matrix [kJ® acts as penalty function 
which constrain the shear-strains Yy^ "to zero as 

plate thickness to span ratio is reduced. Therefore, a 
reduced 2x2 integration rule is normally used for evalua- 
ting shear energy terms to alleviate the overconstraining 
effect, while the full (3x3) integration is adopted for 
nonshear terms. 





CHAPTER 4 


RESULTS AND DISCUSSIONS 

4«1 Preliminary Remarks 

In this chapter, results for stress and stability 
analysis using the present higher order theory for various 
numerical examples drawn from available literature are 
presented and discussedo The reliability of the present 
finite element formulation, programme and applicability and 
necessity of the higher order theory are examined by comparinc 
with the exact 3-dimensional elasticity solutions, closed 
form solutions from higher order theories (denoted as HSDT), 
first order shear deformation theory (denoted as FSDT) and 
classical plate theory (denoted as CPT) results. In the 
various examples, the parameters whose effect is examined 
include boundary conditions, orthotropicity ratio, number of 
layers, side to thickness ratio and fibre orientation. 

The stresses are evaluated at the Gauss points nearest 
to the concerned nodes. Unless stated otherwise the trans- 
verse shear stresses and the geometric stiffness matrix are 
evaluated at the reduced Gauss points. Unless otherwise 
specified 3x3 mesh in quarter plate is used for bending and 
2x2 mesh in quarter plate for buckling. The transverse shear 
stresses are computed using equilibrium equations as well as 
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constitutive relations and presented* In a few cases, the 
effect of different integration schemes is also investigated. 
The computations are done on DEC-1090 system. A note on the 
computational work is presented in Appendix. 

Unless otherwise stated, the deflections and stresses 
are non-dimensionalized as follows : 


w 


wE^h' 


(0,0,0 ) 

^ X* xy' 


(0,0,0 ) h' 

^ X* y* xy^ 

2 




(o , 
^ xz* 


-) h 

yz'^ 




(M^, My, M^y) ^y* ^xy^ p 2 


(N^, Ny, Njjy) ~ ^y* ^xy^ P a 


N 


cr 


cr r; 

E2h 


(4.1) 


The locations for the various quantities presented are 


as follows ; 
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plate at z = 0«5 h 
o’^f cr„ ; Centre of the plate at 2 = 0®5 h 

A y 

^xy* ^xy * Corner of the plate at 2 = 0«5 h 

edge along y-axis at z = 0 
; Mid edge along x-axis at z = 0 

The positive set of stresses and stress resultants is 
shown in Figure 4.1, 


4.2 Presentation of the Results and Discussions 


The results are classified under two ’•rsdi.ngs, as 
follows : 

4.2.1 Stress analysis 

(a) Convergence study, 

(b) Ortho tropic plates, 

(c) Symmetric laminates, 

(d) Unsyrametric laminates, 

(e) Sandwich plates, 

4.2.2 Stability analysis 

(a) Isotropic plates, 

(b) Orthotropic plates, 

(c) Symmetric laminates, 

(d) Unsymmetric laminates. 
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4.2«1 Stress analysis 
4 •2.1*3 Convergence study 

A simply supported 3-layer orthotropic laminates with 
different properties for the middle and outer layers is 
considered with the following data : 

Properties for the middle layer : (in terms of stiffness 

coefficients ) 

C refer a-3i 


^22 

= 0.543103 , 

^11 

^12 

0.23319, 

*^13 

= 0.010776 

fr~ = 0.098276 , 
^11 

^33 _ 

^11 

0.530172, 

^44 

^11 

= 0.262931 

c 

= 0o26681 , 

^11 

^66 

0.159914 




•Oio = = 0.8 h 

12 m 

The properties of outer layer are defined by 


R = = 10 , a = b = 10 in. ^/h = 10, 

^mid 


= 100 psi, = loO 


(4.2) 


The convergence of the non-dimens ionalized maximum 

deflection to the exact 3-D result [44] with mesh refinement 

is shown in Figure 4.2. The convergence of maximum 

o and o are also shown in the same figure. It can be 
xz yz 
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seen that for a 3x3 mesh in quarter plate, the results are 
convincingly accurate and further refinement does not seem 
to be necessary. 

4.2*l.b Orthotropic plates 

A simply supported square orthotropic plate with the 
fo}3~iv:irg properties and subjected to uniform transverse 
pressure is considered for study. 

= 20c83 X 10^ psi, E2 = 10.94 x 10^ psi, 

f% fs 

Gi 2= 6.10 X 10 psi, ^23~ x 10 psi, 

Gi 3= 3.71 X 10^ psi, ^^2= 0*44, 

a = b = 10 in. , = 100 psi (4.3) 

The results are presented in Table 4.1 for a/h ratios 
10 and 20. The comparison is done with Reddy's higher 
order theory results [15j and exact 3-D elasticity results 
[44], The results compare reasonably well. The shear 
stress values computed using equilibrium equations show 
some difference from the values computed using constitutive 
relations . 

4.2.1. c Symmetric laminates 

A 3-layered cross ply (0/90/0) laminate under uniformly 
distributed transverse loading is chosen for study. Thick- 
ness and material of all laminae are same showing the 
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orthotropic characteristics as follows ; 

Ej_ = 19.2 X 10^ psi, E 2 = 1.56 x 10^ psi, 

Gi 2 = =0.82 X 10^ psi, G 23 = 0.523 x 10^ psi, 

■Di 2 = 0o 24, a = b = 10 in. (4.4) 

The edges y = 0,b are considered invariably simply 
supported while x = 0,a are considered SS, CC, FF, SC and 
FS (S, F and C denoting a simply supported, a free and a 
clamped edge respectively). 

For comparison sake, the results for SS, CC and FF 
boundary condition cases are shown in Tables 4.2 to 4.4. 

The results from another HSDT [22], FSDT [22] and CPT [22j 
are presented along with for comparison. The present 
results closely compare with the HSDT results. The stress 
values obtained from the higher order theories and FSDT 
exhibit appreciable difference which are more significant 
when a/h ratio decreases. But, w values computed by the two 
approaches show little difference between them. The classi- 
cal plate theory (CPT) results compare only for b/h = 20, 
when shear deformation becomes negligible. 

Figure 4.3 reveals that the transverse shear effects 
are appreciable for moderately thick to thick plates, 
irrespective of the type of boundary conditions. In 
Figures 4.4 and 4.5 the non-dimensional normal stresses are 
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plotted vs side to thickness ratio and compared with 
FSDT and CPT results. It is observed that for increasing 
side to thickness ratio, the FSDT solution converges from 
below, while the higher order theory solution converges 
from above. 

The variation of stresses through the laminate thick- 
ness is displayed in Figures 4*6 to 4.8 for a/h = 10 for 
varying boundary conditions. It can be seen that the trans- 
verse shear stresses predicted by the present theory are 
non-vanishing at the bounding surfaces. All the same, their 
values are relatively insignificant. As is seen from Fig. 
4.8, the shear stresses computed using constitutive rela- 
tions violate the interface continuity. The shear stresses 
computed using equilibrium equations provide a continuous 
variation. 

4.2.1.d Unsymmetric laminates 

The flexural behaviour of a 2-layered simply supported 
angle ply laminate under uniform transverse pressure is 
studied. The properties and thickness for each layer 
remain the same. The following data is used : 

= 40 X 10^ psi, E 2 = 10^ Psi» ^^2 “ ^ 

= 0.6 X 10^ psi, ^15 = 0.25, a = b = 10 in. 

= 100 psi. 

o 
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The nondimensionalized maximum deflection and stress 
resultants obtained from the present theory are compared 
with results from Ren's closed form solutions of HSDT [18] 
for a/h = 10, 100 and classical plate theory [18] for 
a/h = 100, for different fibre orientations. As seen from 
Table 4.5, the present results compare well with the closed 
form solutions of the Ren’ s higher order theory. For thin 
plates, they converge to the classical plate theory results. 

It is seen from the results that the deflection values 
do not vary appreciably with the fibre orientations. But, 
the values of N^, and drastically reduce to a low 
value for 45°. This indicates that the coupling effect 
is considerably less for 00 = 45°. 

4.2. l.e Sandwich plates 

Example 1 : A simply supported square sandwich plate with 
aluminium faces and balsa core subjected to uniform trans- 
verse pressure is studied with the following face properties 

= 10 X 10^ psi, t^ = 0.032 in. (4.5) ' 

The objective of the study is to establish the applica- 
bility of the present theory to the analysis of sandwich 
plates, by comparing with available test [45] results. The 
results are presented in Tables 4o6 and 4.7 for varying 
spans, core thickness, core shear rigidity and load inten- 
sities. The present results compare satisfactorily with 
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the test results. Since all the plates considered therein 
fall within the thin plate domain, only slight improvement 
is found in some of the cases over the first order theory 
results [ 45 ] (also presented). It is further observed that 
when the core shear rigidity is low, the consideration of 
face shear rigidity yields better deflection values in 
majority of the cases. Though this does not lead to a clear 
conclusion in this regard, because of the number of parameterj 
varied, it prompts us for a more systematic investigation, 
which is done in the next example. 

Example 2 : A simply supported square sandwich plate with 
the following properties is taken for study : 

^ ^ 13 

Face : = 25, Eg = lo'" psi, 


G 

= 0o2, = 0o25 . 


G 


23, 


Isotropic core : E = 0, G is defined by R — ^ 


a = b = 10 in., = 100 psi, a/h = 10 


(4.6) 


mesh ; 2 X 2 in quarter plate. 


The effect of considering the shear rigidity of faces, 
on the computed response is studied for varying core to face 
thickness ratio (C/t) and ’R’ values. The results are 
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displayed in Table 4.8o It can be noticed that as c/t 
ratio increases, for any ’R’ value, the effect of considering 
shear rigidity of faces diminishes. This is because, as the 
face becomes relatively thinner, the transverse shear effects 
in the face become less and less significant. Secondly, it 
is predominantly seen that as R increases, i.e., as the 
transverse shear rigidity of core relative to the shear 
rigidity of face decreases, the consideration of face shear 
rigidity makes a marked change in the results. This is more 
so for thicker faces. The explanation can be that when the 
core is relatively weak in shear, the resistance of the faces 
to transverse shear deformation become relatively more signi- 
ficant and effective. The effect of considering the shear 
rigidity of faces is not predominant when core shear rigidity 
is high, as due to the fact that the faces are relatively 
thin and their shear resistance will become insignificant 
compared to that of the thick, rigid core. 

Example 3 : A simply supported sandwich panel with thin 
faces is taken for study with the properties mentioned 
below, 

face ; = 10 X 10^ psi, -Di. = 0.3, t^ = 0.005 in. 

core I = 10,000 psi, t^ = 0.4 in. 

c ^ 

^ 1.0 xn « 


(4o7) 
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The maximum deflection and stresses are computed using 
the present theory and shown in Table 4»9 along with results 
from Allen [46], The results compare well for both aspect 
ratios* The effect of integration scheme on the computed 
values is also studied. It is found from the Table 4.8 that 
for the relatively thick plate that is considered in the 
problemo, the order of integration does not have a pronounced 
effect on the values, for the element under consideration. 

Example 4 : Rectangular sandwich plates with laminated faces 
of size le in- are considered for study. The data 
used follows : 

CFRP face plate : 3 layers each of thickness o-oos'in. for 

each face plate. 

Gi 2 =o t<i.-.cfpsr, = 0»327 

aluminium honey-comb core : 

Gi 3 G23 = (iii , 

p =s o • I U- £.2. psi 

o 

Four types of lay ups are considered for the specimens 
for which experimental results are available. 

SPl ; [303/core/303], SP2 : [03/core/03] 

SP3 : [30/-30/30/core/30/-30/30] 


SP4 : [0/90/0/core/0/90/0j 
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The core thickness for SPl and SP2 is 10 mm, while 
for SP3 and SP4, it is 7 mm. (4.8) 

The results are displayed along with the experimental 
L47] and FSDT results [47] in Table 4.10. The deflections 
are computed for both simply supported (SS) and clamped (CC) 
boundary conditions, as the test boundary conditions are 
stated to be in between the two types, but intended to be 
clamped. It is seen that the higher order theory results 
for the clamped condition are closer to the experimental 
results than the FSDT results. But, there is considerable 
difference between the experimental results and the present 
results. The reason is that the clarnpsd boundary conditions 
were not satisfactorily achieved during experimentation, as 
stated in the reference. [ 47] . 

4.2.2 Stability analysis 
4. 2. 2. a Isotropic plates 

A simply supported, square isotropic plate is consi- 
dered with properties and other data as given below : 

E = 2,5 X 10^ psi, ‘D= 0.25, a = b = 10 in. (4.9) 

The uniaxial buckling loads are computed. The results 
are displayed in nondimensionalized form in Figure 4.9 along 
with exact 3-D elasticity [35] and classical plate theory 
[35] results. The figure exhibits the close agreement of 
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of the results computed using the higher order theory with 
the 3-D elasticity results for range of plates from thick 
to thin. The mesh used for computing the values being 2x2 
in quarter plate, it can be expected that as mesh is refi- 
ned, the agreement with the exact solutions will be more 
remarkable. It is clearly seen that as the plate thickness 
increases the buckling parameter decreases. This is because 
of the transverse shear deformation that becomes predominant 
for thicker plates. The classical plate theory which neglect 
the transverse shear effects is found to be erroneous for 
thicker plates . 

4.2.2.b Orthotropic plates 

A simply supported square plate with the following 
orthotropic properties in the form of ratios of stiffness 
coefficients is treated for study : 


f22 

^11 

= 0 .543103, 

^12 

= 0.23319, 

^13 

= 0.010776 


^11 


^11 


^23 

^11 

.= 0.098276, 

^33 

^11 

= 0.530172, 

^44 

^11 

= 0.262931, 

^55 

Cll 

= 0.26681, 






a = b = 10 in . , 


*^11 - 


(4.10) 
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The nondimensionalized uniaxial critical buckling load 
is plotted versus a/h ratio in Figure 4.9. The exact elasti- 
city [44], Mindlin’s [44] and CPT [44] results are also shown 
for comparison. Once again, N decreases with increase in 

Cr. i 

a/h ratio which can be attributed to the influence of the 
transverse shear deformation. The thin plate theory is seen 
to provide optimistic values of buckling load, the errors 
increasing with thickness. The thin plate assumptions 
increase the stiffness of the structure, thus yielding higher 
values . 


4.2.2.C Symmetric laminates 

A simply supported symmetric cross ply laminate is 
studied. The following data is used : 


'13 


'12 


0 . 6 , 


G 


23 


= 0.5, '0 


12 


0.25, = 10^ psi. 


a = b = 10 in., a/h 


10 


(4.11) 


The nondimensionalized critical buckling parameter is 
evaluated for different number of layers and different 
degrees of orthotropy. The values are shown in Table 4.11. 
The comparison is done with a HSDT [39], FSDT [39] and 
CPT [39] results. The present results match closely with 
the HSDT result. For the moderately thick plate that is 
studied, the response is not significantly different from the 
FSDT results. Probably, the effects of the assumptions of 
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higher order theory can be felt for thicker plates. It is 
observed that the CPT results overestimate the buckling load. 
The error increases with increase in ortho tropicity ratio. 

It can be noticed that the integration rule adopted has 
very slight effect on the results. Still, the results 
obtained from full integration are consistently slightly 
higher than their counterparts. This is due to the over- 
stiffening effect when the shear terms of the stiffness ma- 
trix are formed using full integration. 

Figure 4.10 displays the variation of buckling para- 
meter with degree of orthotropy. The buckling parameter is 
found to increase with the degree of orthotropy. Further, 
from Table 4.11, it is seen that the number of layers does 
not appreciably affect the critical buckling load. However, 
for higher ortho tropicity ratio, there is a slight noticeable 
increase in the values with increase in number of layers. 

4, 2, 2. 3 Unsymmetric laminates 

Example 1 : A simply supported square antisymmetric cross 
ply laminate is considered for study, with the properties 
of individual laminae and other data as in 4.11. 

The Table 4.12 is devoted to the comparison of the uni- 
axial buckling parameter (N ) obtained from the present 
approach with those available in literature, and to emphasise 
the effect of inplane orthotropicity ratio and number of 
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layers® The present results compare better than the 
Reddy's HSDT results [37], with the exact 3-D elasticity 
results by Noor [48], especially for higher degrees of 
orthotropy. The numerical results allow one to infer that 
for a thick plate, and when material exhibits high ortho- 
tropicity ratio, the effect of transverse shear deformation 
needs to be incorporated* It can also be remarked that the 
difference betvi^een the classical buckling load on one hand 
and their shear deformation theory counterparts on the other, 
increases drastically for high orthotropicity ratio, as the 
number of layers increase. The values, in general, 
increase remarkably with increase in number of layers from 
2 to 4, but thereafter, marked increase is not noticed. This 
behaviour is due to the fact that, when there are two layers, 
the coupling between extension and bending is maximum, 
leading to a decrease in the buckling load. With increase in 
number of layers to four, the coupling effect is significant!'^ 
reduced, thereby increasing the buckling load. With further 
increase in number of layers this reduction is relatively 
less drastic. A comparison between Tables 4.11 and 4.12 also 
reveals the role of coupling in bringing down the buckling 
strength, the effect becoming less significant for larger 


number of layers 
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Example 2 ; A simply supported, square, antisymmetric angle 
ply laminate with the following properties is studied for 
response to uniaxial uniform compresive loading. 

= 10^ psi, =0*6, G23/E2 = G31/E2 = 

^12 “ 0*25, a = b = 10 in. , Numbers of layers = 2 (4.12) 

The uniaxial buckling parameter N for varying fibre 

orientation and degree of orthotropy are computed and shown 

in Tables 4<.13 and 4.14 for a/h = 10 and 100 respectively. 

The results match well with the results from available 

HSDT [41] for all fibre orientations and orthotropicity 

ratios. The values increase with increase in inplane 
cr 

orthotropicity ratio for all fibre orientations, and side to 
thickness ratios. For © = 5°, this trend is most pronounced, 
as the fibres being closer to the axis of loading, any 
increase in the strength along the fibres will be more signi- 
ficantly reflected along the loading axis. A comparative 
study between the Tables 4.12 and 4.13 reveals that the 
buckling parameter values decrease with increase in thickness 
the shear deformation effects being the cause. 



CHAPTER 5 


CONCLUSIONS 


5*1 General 

A higher order shear deformation theory for the bending 
and buckling analysis of laminated composite plates is 
presented. Certain general conclusions can be drawn from its 
application to numerous examples. 

* The displacement based finite element formulation for the 
buckling analysis of laminated ccniposite plates using a higher 
order displacement model had not been attempted so far in the 
literature. The various numerical examples worked out indicate 
the reliability of the present formulation. The close agreement 
of the results with exact solutions makesit acceptable for 
general use in the analysis of the stability of laminated 
composite plates. 

The formulation for the stress analysis has been tested 
for various boundary conditions in Section 4.2.1i.Co Also, 
in Section 4,2oi«d, examples worked out for unsymmetric 
laminates for varying fibre orientations and side to thickness 
ratios have close agreement with closed form solutions of higher 
order theory. The consistently accurate results lead to the 
conclusion that the present theory is applicable to a wide 
range of problems associated with bending of laminated composite 
plates . 
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* The adaptability of the present formulation to sandwich 
plate problems is thoroughly investigated in Section 4*2. l.e. 

The parameters involved include varying core and face properties 
and various combinations of fibre orientations for the CFRP 
faces. The satisfactory results obtained allow us to be confi- 
dent about the predictive capability of the higher order theory 
with regards to stress analysis of sandwich plates. 

* In Section 4.2.1,e, example 2 was devoted to a systematic 
examination of the effect of considering the shear rigidity of 
faces in sandwich plate analysis. The study reveals the 
follows : 

(a) For sandwich plates with thin faces (c/t > 50), considera- 
tion of the shear rigidity of faces has no significant effect, 
even with relatively weak cores. 

(b) For thicker faces, when the shear rigidity of faces is high 
relative to that of the core (R > 10 approximately), the neglect 
of face transverse shear leads to high, erroneous deflection 
values . 

(c) The moment resultant values are generally not much affected 
by the face shear considerations. 

In general, it can be concluded that the face shear rigidity 
is to be accounted in the deflection analysis of sandwich plates, 
except for the case of very thin face and cores with high shear 
rigidity. 
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* As can be seen from the results in Tables 4o2 to 4.4 , 
the first order shear deformation theory is satisfactory for 
the defloctioii calculations. But, for thick laminates, the 
stress values from FSDT are highly erroneous and it becomes 
inevitable to resort to the higher order shear deformation 
theory. The classical laminate theory can be relied upon in 
the thin plate (a/h > 20) domain only. Even there, it lacks 
the capacity to predict the transverse shear stresses. 

The FSDT results for the uniaxial critical buckling load 
for symmetric cross ply laminated plates is given in Table 4.11, 
For the type of laminate studied, it can be seen that the first 
order shear deformation theory will suffice to determine the 
lowest buckling load. Again, the classical plate theory is not 
reliable for thick plates, the error increasing with degree of 
orthotropy , 

* It is clear from Tables 4,9 and 4,11 that the type of 
integration adopted has little effect on the results when the 
9-noded Lagrangian element is used either for the stress analysis 
of sandwich plates or buckling response of a laminate. 

5.2 Reco'ir.cndaticns for Future Work 

Some of the future work that can be carried out as an 
extension or improvement over the present investigation are 
suggested below. 
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some effort, this work can be extended to the stress 
analysis of stiffened laminated plates and laminated composite 
shells. The response to combined transverse and inplane loading 
has not been studied in the present investigation. This can be 
included in a future work. The present work can form the basis 
for a formulation for the nonlinear analysis. 

The buckling formulation and program that is developed can 
be utilized to analyze the stability of sandwich plates. The 
program also caters for biaxial and shear buckling of laminated 
composite plates. These could be investigated in a future work. 
Formulation can be developed for buckling due to combined trans- 
verse and inplane loads including shear. 

An improvement over the buckling formulation presented is 
to consider all the components of the nonlinear strains such as, 

^ du 1 , /dv%2 . /dw\2-j ^ 

®x = ^ 2 + (3^} J, etc., 

in forming the initial stress matrix. The problem becomes more 
involved, A stress analysis may have to be carried out to 
determine the stress distribution and the results used to form 
the initial stress matrix. 
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Table 4.3: Maximum deflection and stress in a 3-layered crossply 
laminate under uniform transverse pressure. 
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■'■'TATION 


A FORTRAN impl^rrietitation of the present formulation is 
used to solve a variety of r,.^;nsrical examples. Programmes have 
been separately written for bending and buckling cases. Sepa- 
rate routines are written to generate the stiffness matrices 
for symmetric and as. y I'. tric laminates and stored in separate 
files, any of which may be run along with the main file depend- 
ing on the problem. The separate formulation and program for 
the symmetric laminate rings down the CPU time to a quarter of 
that for an uasyr-r-chric laminate. 

In view of the T'r'j : matrices that would be involved, the 
global, stiffness matrices are stored in a compact form using 
the profile (skyline) storage technique. In this technique, 
coefficients of the global stiffness matrix within the non- 
zero profile of the matrix only are stored. This always requires 
less storage than a banded storage. Secondly, the storage re- 
quirements are not severely affected by a few very long columns. 

The procedure involves the calculation of the number of i 
equations to be solved after the equations with | 

prescribed displacements. The equation number for each degree ; 
of freedom of each node is assigned. Then, using the elemental 
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connectivity and equation numbers, the ccl'jnu heights above the 
principal diagonal are calculated. Using these column heights, 
positions of diagonals in the single subscript array for the 
stiffness matrix storage are determined. This acts as a painter 
array with the aid of which the actual location of any other 
coefficients stored in the compact array can be decided. 

The Gauss elimination is used for the solution of the 
simultaneous equations that arise in the stress analysis* The 
subroutine COLSOL available in Ref. 43 which operate on the 
stiffness matrix stored in compact forms is used. 

The solution for the eigenvalue problem that arises in 
the stability analysis is obtained using the subspace iteration 
method. The subroutine SSPACE listed in Ref. 43 is used to 
implement this method in conjunction with the profile storage. 

The routine yields the specified number of lowest eigenvalues 
and the corresponding eigenvectors 
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